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An increase in the size of coherent domains in the one component field theory under the 
influence of a uniformly changing external magnetic field near the critical end-point T<i, = Tc, ftc> =0 
was proposed recently as an estimate also for the variation of the chiral correlation length of QCD 
near its respective hypothetical end point in the Tqcd — Hqcd plane. The present detailed numerical 
investigation of the effective model suggests that passing by the critical QCD end point with realistic 
rate of temperature change will trigger large amplitude oscillations in the temporal variation of the 
chiral correlation length. A simple mechanism for producing this phenomenon is suggested. 



I. INTRODUCTION 



o 
o 

^, The quality of our knowledge of the phase structure of QCD at high temperature and finite baryon density is an 
• important benchmark for our understanding of strong interactions. 
m , A critical end-point of the first order phase transition hne in the T — /x-projection of the QCD phase diagram was 
■ conjectured Q to follow from the compatibility of the following observations: 

' i) Lattice simulations |^ indicate that the phase transformation at zero chemical potential with realistic quark 
masses is a crossover, characterized by an analytic variation of the thermodynamical potential with the temperature; 
Cn| ii) At zero temperature there is a first order phase transition from the hadron phase to more exotic phases as a 

^ function of fi j^] which continues as a first order line into the (T, /i)-plane. 

2^ ' For the strong matter the significance of this end-point would be similar to that of the Curie point for ferromagnets. 

Recently, some important progress was realized in the search for the {^EiTe) location of the end-point with non- 
perturbative lattice studies There is a chance in current heavy ion collision experiments that it can be observed 
\ experimentally, since with increasing collision energy/nucleon the central rapidity particle spectra explore regions of 
T-H 1 the phase diagram corresponding to decreasing chemical potential. An obvious class of the signatures would reflect 
the increasing size of the coherent fluctuations in the order parameter cr-field in the neighborhood of the end-point 
|§| • In this region mostly the a-field will be excited, since its mass is the lightest near the critical end-point in view of 
O i' the amount of explicit chiral symmetry breaking which keeps the pions massive. This coherence should be reflected 
I by the statistics of the main decay products of the cr-field, the pions. 

Starting from equilibrium at some Tq > Tc, and passing with finite velocity near the end point, the system unavoid- 
(-H ably slows out of equilibrium. In contrast to the equilibrium characterization of the second order transition, finite 
• • . maximal correlation length is reached with a certain amount of supercooling. A substantial increase in the correlation 
' length will be a reliable signal for the existence of the critical end-point. 

. Quantitative predictions for this phenomenon should rely on non-equilibrium field theory. For the moment it is 
' hopeless to simulate directly the far from equilibrium behavior of QCD matter. One recognizes however, that the 
5^ , mass of the a in this region is separated by a gap from the other mesonic excitations. Therefore, one arrives at 
the conclusion that an effective dynamical theory of the longest wavelength excitations in this region is in the same 
universality class as the Ising model. At present no quantitative matching is known between the original theory and 
its one-component effective version. Still, universal features of the class A of critical dynamics (in the classification 
of Hohenberg and Halperin [^) are expected to occur. ^ 

Recently Berdnikov and Rajagopal (BR) proposed an intuitive mapping. They approximately identify the 
magnetic field {h) of the Ising model with the temperature of the QCD (the Ising reduced temperature (r) axis is 
nearly parallel to the chemical potential axis of QCD). They have checked that the results are not sensitive to a 
moderate tilt in this mapping. Next, they have proposed a dynamical equation which describes the evolution of the 
inverse correlation length (the mass of the a meson) when the system passes through the critical end point of the 
Ising model with finite velocity under different angles. This equation was shown to depend on a single non-universal 
parameter, proportional to the rate of change of the Ising magnetic field. 



*e-mail: mazsx@cleopatra.elte.hu 
^e-mail: patkos@ludens.elte.hu 
■'■e-mail: sexty@cleopatra.elte.hu 
^e-mail: szepzs@cleopatra.elte.hu 

A more complete theory reflecting the direct infiuence of the temporal variation of a conserved baryon number density on 
the a field falls into class C |7|, and will be subject of future study. 
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Finally, a semi-quantitative correspondence was proposed between the relevant QCD temperature range (Tq = 
180 MeV, Tfreezeout = 120 MeV) and the dimensionless strength of the Ising magnetic field (/iq = —0.2,hfinai = 0.1). 
The non-universal constant was varied in a wide range, since it relates essentially the rate of cooling of the QCD-matter 
to the speed of the variation of the external field h. The main result of Rcf. |j] was a prediction for the variation of 
the order parameter (e.g. a) correlation length with the variation of h in the neighborhood of the critical end-point. 

The variation of the correlation length during the passage through the critical point was compared to its value 
taken at the non-critical starting value, which corresponds to the equilibrium energy density of the QCD plasma 
produced in the collision {T = Tq). It is the relative increase in the correlation length at the freeze-out of the system 
(T = Tfreezeout), where most of the pions are coming from, which is the most important issue. A factor of 2-3 increase 
was signalled, which is claimed to be observable under the accuracy of the present heavy ion experiments. 

The aim of our present investigation is to check the accuracy of the first order relaxational dynamics assumed for 
the inverse correlation length. For this we have integrated exactly the equations of motion of the one-component 
classical three-dimensional scalar field, <&ti(x, t) on lattices of size TV = 16 — 128: 

$d(x, t) - (A - m^) $rf(x, t) - ^'^di^, t) - hd{t), (1) 

and measured simultaneously the variation of the order parameter and of the correlation length. 

Though the hadronic system freezes out at Tfreezeout ~ 120 MeV which corresponds to h — 0.1 in the effective 
system, we have followed the variation of the correlation length and of the order parameter to higher values of h. 
This enabled us to recognize the relevance of a second order dynamics in the effective equation of motion of these 
characteristics. 

It turns out that the order parameter (OP) obeys an equation which is slightly more complicated than the one 
proposed by Halperin and Hohenberg for this class. It is formally analogous to the differential equation of a damped 
oscillator. In order to achieve good quantitative description of the OP-trajectory obtained from the numerical simu- 
lations, one has to take into account the effect of slowing out from equilibrium while the system passes by the critical 
point. 

Our paper is organized the following way. In Section II the method of the numerical solution of is shortly 
outlined. The methods of analyzing the time evolution of the system are presented in more detail with special 
emphasis on the determination of the spatial correlation length. In Section III the trajectories of the most important 
quantities characterizing the state of the system are mapped out in the (r, h)- and the (OP, /i)-planes when passing 
at different distances by the critical Ising end point. In Section IV we present the results for the non-equilibrium 
/i-evolution of the correlation length for several values of = dh/dt and compare them with the estimates of BR. 
Finally, in Section V a quantitative phenomenological interpretation of the measured order parameter trajectory is 
offered. Our dynamical description is compared to the equation of motion proposed for this class by Halperin and 
Hohenberg. Based on the proposed effective relaxational dynamics we suggest also a simple way to account for the 
variation of the correlation length. Conclusions of our investigation are summarized in Section VI. 



II. METHODS OF SOLVING AND ANALYZING THE EXACT CLASSICAL EQUATIONS 

For setting up and solving numerically the theory we used techniques similar to those applied in our previous 
paper Here we shortly outline the procedure. 

As a first step one has to rewrite Eq.([^) using dimensionless quantities (not having index d), defined as follows: 

t^td/ax, x^Xd/a^, 
$ = y^a^$d, h = yf^alhd, (2) 

where ax stands for the lattice spacing, the powers of which we used to scale dimensionfuU quantities. The dimen- 
sionless mass parameter of the theory |m| = aa;|md| was set to unity. The field equation in dimensionless variables is 
of the following form: 

<^n{t + at) + $„(t - at) - 2<i>„(i) - a'Jal 5](<i>„+.(0 + $„_-(0 - 2<f„(i)) + a?(-<l>n + - h{t)) = 0. (3) 

i 

There were a number of simple quantities routinely monitored in each run. The first was twice the average kinetic 
energy per site (called kinetic temperature also for non-equilibrium field configurations): 

Tk^n{t) = ^J2^l{t), (4) 
n 

the second, the trajectory of the homogeneous (OP) mode: 
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n 



(5) 



Also its fluctuation 




(6) 



was used. 



The thermalization algorithm, which led to the initial state, consisted of two steps. First we have set the initial 
energy density by continuously comparing the desired and the measured kinetic temperatures. Depending on the 
deviation from the targeted temperature an artificial friction or anti- friction term was introduced into Eq.(H). After 
having reached the required kinetic energy density, in a second step the original form of the discretized nonlinear field 
equation was iterated until thermalization was complete. Both steps were repeated until the final kinetic temperature 
was just what we desired. 

The critical temperature T^^c was determined at /i = by locating on the axis the maximum of S^^, or of the 
specific heat, and by finding the point separating zero and non-zero expectation value regions of OP, $. The reduced 
temperature r = {T^ — T^^c)/T^,c was measured relative to T$_c = 1-5 for L = 32 and T$^c — 1-57 for L — 64. 



The main goal of this Paper is to give a quantitative interpretation to the dynamical behavior of the correlation 
length as the system passes by the critical end point. Therefore we need an accurate measurement method for this 
quantity, which is reliable in a dynamical system too. 

1. The most straightforward way is to use the definition of the correlation function: 



Here the overline refers to the spatial average of some quantity at fixed i in a single sample, while the 
brackets stand for the average over the initial conditions. After checking that C(A, t) truly behaves as 
~ cosh[(A — L/2)/^i], one can extract the correlation length ^i. We refer to this characteristic length as the 
"direct" length of correlation below. 

2. Average linear size of the regions of same sign deviations from the space average $ can be taken to estimate the 
characteristic size of coherent fluctuations. Consider $(x, t) — At a given time a histogram was constructed 
by scanning through the lattice for the number of occurrences of site sequences with the same sign deviation 
from $(t), and of a given length A, parallel to the lattice axes. The histogram showed perfect exponential 
dependence on A. Its characteristic length defines the "grain size" ^2- Repeating this measurement for every 
configuration during the time evolution one obtains the function ^2{t). 

The correlation lengths defined by the above algorithms arc different, of course. One expects, however, that both 
definitions capture the same feature of a field configuration and there exist simple functional relationships between 
them. In order to find the relation of to the standard correlation length (or its inverse, the mass) , which is usually 
measured with method ^ (see below) , we studied the equilibrium systems at different values of the reduced temperature 
r e (0, 0.27). Using the algorithms described above one finds the functions ^i(r), as well the standardly used meff{r). 
The elimination of r leads to the relations rrie// = gi{£,i), well approximated by second order polynomials. In principle, 
one should establish this relation for each value of h separately, but in practice the relation found for h = Q gave good 
agreement between the results of the different methods, also away from equilibrium. 



An algorithm for the reconstruction of the effective potential from the real time evolution of a scalar field was 
presented in An effective equation of motion was fitted to the temporal variation of the order parameter. Its 
"force" term was interpreted as the derivative of the effective potential with respect to the field. Then it was easy to 
identify the effective squared mass. 

This time we improved further this algorithm. All spatial Fourier modes of the system are used for the reconstruction 
of the dispersion relations. Three stages of the implementation were worked out. 



A. Determination of the correlation length 



C(A, t) = mx, y,z + A, i)$(z, y, z, t)) ~ (mf. 



(7) 



B. Spectral determination of the mass 
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3. We assumed the validity of the following approximate equation for short time intervals: 



A$ + m;:^^.$ + eA^$ = 0, (8) 

where A is understood as lattice discretization to the Laplacian. Comparing this to the real nonlinear evolution 
we could fit ?7i^jy,e for the short intervals in question in the following way. A spatial FFT algorithm was 

applied to the field configurations generated by Eq.(^. The temporal trajectory of each k mode was fitted to 
the equation dictated by the Fourier transform of Eq.(|8|). The coefficients of the polynomial of k'^ determine 
also the time dependent effective mass. (Here and in the following ki — 2(sinfci/2) stands for the dimensionless 
lattice momentum.) A lower bound for the length of time interval in which Eq. is fitted comes from the 
consistency criterion that the time interval we averaged over must be much greater than the resulting l/nieff 
time scale. For out-of-equilibrium field configurations there is also an upper bound for the time interval which 
is set by the variation rate of the parameters r, h. 

4. For the two-point function related to the OP-susceptibility j|] the relation 

T 



= Z-'{miff + k' + ek^ + Oik'')) (9) 



holds in equilibrium llOl . Replacing the ensemble average in the denominator of the left hand side by averaging 
over the k- modes characterized by equal fc^ and using T^.^^^ ^,2 — |*i>kp for T we found that our system obeys 

Eq.(^ for k^ > 1 without any time averaging. With these replacements we could fit the value of ■m'^ji'{h) 
continuously both near and far from equilibrium. The wave function renormalization constant was found to be 
equal to 1 up to 0.5 percent in all cases. The coefficient e of the fourth derivative correction term was checked 
to be neghgible. 

5. Since most of modes with fc^ > 1 nicely follow the above behavior, one might improve the statistical confidence 
of the above algorithm by summing over the contribution of all modes (including k'^ < 1), which yields the 
approximate equality 

with denoting the field fluctuation as defined in Eq.(^). Its measurement does not require the application 
of the time consuming FFT. The value of the function on the right hand side can be tabulated and Eq. (p^) is 
quickly solved with help of a look-up table plus interpolation for every measured value of 5^^. 

In addition, as a quick reference the so-called Hartree squared mass estimate wl^artree = ^1 + 3($^ + 6^^) was 
also used. It proved useful in interpreting qualitatively effects in the motion of modes with different k and apparently 
related to the presence of an /i-dependent common effective mass. 

After careful testing of the simplified procedures against the conceptually better funded methods in equilibrium we 
decided to use Algorithm || throughout this paper. 

Although the different algorithms were normalized to yield equal masses in equilibrium, it is not obvious that 
they will agree also for a non-equilibrium crossover connecting two points of the (r, h) plane. We shall return to the 
comparison of the non-equilibrium results obtained with different methods at the end of section IV. 



III. PASSING BY THE CRITICAL POINT 



We have solved Eq. (||) , the discretized version of Eq. ([I]) numerically. A starting configuration was thermalized in 
presence of the initial magnetic field {ho = —0.2), in such a way that a predeflned value of the average kinetic energy 
density was reached. This thermal state was taken for the initial configuration when solving the field equation with 
an external magnetic field tuned with constant velocity dh/dt. With the parameters we used the correlation length 
of the initial configuration was approximately equal to one lattice unit. 

In the present investigation the same range of the parameters r, dh/dt was covered as in [Q. For each value of 
r,L,a= (dh/dt)^^ runs with ~ 20 different equilibrium configurations were averaged. 

An approximate idea about the part of the Ising phase diagram explored in our numerical investigation can be 
given by drawing the measured (r, h) trajectories for different values of r initial and a (see FigJ^). 



4 




FIG. 1. The measured trajectory in the (r, /i)-plane for two different initial temperatures (sohd hues). The crosses show 
the reduced equihbrium temperature r{h) corresponding to energy density e{h). The size of the system was L = 64, and the 
inverse rate of change of h: a = 100. 



Due to the time dependence of the external force h{t) the fuU energy density e{t) is not conserved, we shall 
parametrize it as e(/i). In principle, it might depend on the parameter a too, but for its range investigated in this 
paper, we did not experience any a-dependence of the energy density. The reduced (Ising) temperature, which is 
defined through the average kinetic energy per site, slowly drops in the first part of the crossover. However, having 
passed by the critical point {h — 0, r = 0) a, deterministic oscillation in r starts, which survives the averaging 
over the ensemble of initial thermal configurations. For comparison we have displayed in Fig.|l| also the reduced 
temperatures corresponding to the equilibrium states belonging to the different values of the energy density e(^). The 
latter was determined by stopping the variation of /i at a certain value hstop, and then relaxing the system with fixed 
h = hstop (and hence, with conserved energy) to equilibrium. After the thermalization of this state was complete, we 
have measured its equilibrium temperature. These runs were performed for every non-equilibrium trajectory for 40 
equidistant intermediate values of hstop G (—0.2,0.2). 

Fig.|l| clearly shows how the system slows out of equilibrium. As the critical point is approached, the thermalization 
time scale grows above the time scale ('^ ha) of moving in the phase diagram. Thus, in the second part of the crossover 
two effects seem to be present: First, the spectral density of the configuration will correspond to the equilibrium at an 
earlier h value - this means an overcooling in QCD language. Second, the low-k modes get highly excited and begin 
to oscillate with a frequency determined by an effective mass scale. The oscillation is synchronously present in the 
entire spectrum, i.e. the kinetic energy of each "well-behaving" (fc^ > 1) mode oscillates coherently. The oscillation is 
not due to any direct strong coupling between the modes, but is driven by the oscillation of the homogeneous mode 
($) (see Fig p|a). The UV modes adiabatically follow the slow oscillation of the order parameter(see Fig.^)), which 
enters parametrically through a common effective mass in the respective equation of each mode. 




a) b) 



FIG. 2. a) OP trajectories for different dh/dt = a ^ values (L = 64, riniuai = 0.15). The equilibrium ^(/i) curve is 
represented by crosses, b) Temporal variation of the left hand side of (^) synchronized by the OP oscillation. 

One may define the error of OP trajectories from different thermal initial conditions as the standard dispersion of 
the OP(t) values at fixed t. This error comes out at « 0.007 for /i < and w 0.02 for h > 0, which is about 100 times 
smaller than the amplitude of the oscillations. 
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For comparison the equilibrium OP values are also displayed in Fig. ^a. These values come from thermal solutions 
of Eq.(jl|) with the (r, h) values chosen from the equilibrium points of Fig.|l|. 

The OP-evolution displayed in Fig.||a is by itself a challenge seeking quantitative interpretation. We shall elaborate 
on it in detail in Section 4. 



IV. HOW LARGE THE CORRELATION LENGTH GROWS? 



The main physical motivation of our investigation was to answer the question in the title of this section. We have 
used the fast method described in Section 2 for deducing the actual non-equilibrium value of the correlation length, 
that is the inverse cr-mass. 




FIG. 3. Evolution of correlation length during a crossover compared to the equilibrium values along the (r, h) route depicted 
in Figjl|. (L = 64, r^.^) — 0.083, ri^i,) = 0.15). Crosses give the results of equilibrium measurements. 

In Fig.^a we display the /i-history of the correlation length obtained on a 64'^ lattice with an initial reduced (Ising) 
temperature = 0.083 for four values of the a parameter. The (r, fi) trajectory starting from this passes the closest 
by the critical end point (r^m — 0.007), see the lower curve in FigJ^. Again, the equilibrium correlation length values 
obtained from analyzing long thermal time evolutions with an identical method, are displayed in the same plot. 

When compared with the estimate of the most dramatic difference is obvious: even the slowest variation of h 
considered in BR is in reality much too vehement. The critical slowing down of the internal interactions pushes the 
system far out of the equilibrium state corresponding to the actual magnetic field value. This is the reason that after 
passing the maximum of the correlation length, reached with the expected "supercooling" in h, ^ starts to oscillate 
and the correlation length very steeply drops to a minimum. This oscillation has opposite "phase" if compared to 
the (Ising) temperature. A probable explanation to this phenomenon is, that a shrinking correlation length means 
larger cr-mass and hence an increase in the frequency of the microscopic oscillations of the UV-modes. The energy of 
a weakly coupled UV mode — like of a tuned linear oscillator — gets larger with the increase of its frequency. 

The qualitative picture is the same in a wider neighborhood of the critical end-point as one can see from Fig.||b 
which corresponds to a trajectory passing somewhat farther (r^ = 0.15, the upper curve in Figjl]). 

The period of the oscillation in QCD temperature based on the Tqcd — h correspondence proposed in |7| seems 
to be much smaller than the spread of the freeze-out temperature estimates appearing in the literature. Therefore 
our main qualitative result is that the expected increase in the coherence of the pion radiation might be missed if 
the actual freeze-out would happen at a temperature slightly beyond the maximum of the correlation length. Ideally, 
accurate measurements of the freeze-out temperatures of different meson species coupled to a might allow to map out 
the variation of its correlation length across the hypersurfaces of the respective "last scatterings" as predicted by our 
analysis. 

One can compare the exact non-equilibrium ^{h, a) function to the estimate of BR in some more quantitative detail. 
The maximal values of the correlation length in units of the initial are in the range (2.5 — 3.5)^o- The amount of 
supercooling is generically larger from the numerical solution of the <I>^ dynamics, in comparison to the result of the 
first order dynamics conjectured in 0. 

A good measure of the amount of the physical supercooling as a function of the /i-velocity is offered by the shift in 
the location of the first maximum of the correlation length, h{^,nax)- This value turns out to scale with the velocity 
of the /i- variation: 

MUa.,a)=m-°-™i. 



6 



The values of £,max themselves follow also a scaling form as suggested in ||7| on the basis of dynamical scaling consid- 
erations. We find numerically 

e™a.(a) = c'a°-2"±0-0i. (12) 

The value of this exponent agrees with the prediction of BR, when it is applied to the class A. Then their prediction 
for the exponent yields v/(35/(l + zv/135) = 0.222, with = 0.630, (3 = 0.326, 5 = 4.8 and z = 2 + (6 ln(4/3) - 1)tj, 77 = 
0.0335 When the initial (Ising) temperature is tuned to approach the critical point the closest by 2-3% from 
above, it turns out that the form of the ^(/i) function qualitatively remains the same as described above. The ^ values 
at the first maxima are about four times larger than those taken in the first minima following them. 

The second maxima of ^ appearing in our numerical solution very probably cannot be observed experimentally, 
since the fireball breaks up into non-interacting mesons before reaching the corresponding low temperature. 

The sensitivity of the results to the size of the system is also an important issue, since at present the linear size of the 
plasma droplet is estimated to about 6fo- Since we have chosen ^0 for the lattice constant, one might expect to reach 
the maximal allowed correlation length in a much smaller lattice volume already. We have tested the robustness of 
the above conclusions by varying the lattice size between 16 and 64. No important variation was seen when changing 
the size from L = 64 down to L = 32, but a 20% drop in the maxima of ^ appears when going down to L = 16. 

The lattice we used in this investigation is rather coarse (|m|aa; = 1). The idea behind this choice is that we work 
near the critical point, in the scaling regime. At equilibrium when the correlation length grows very large it is obvious 
that the actual value of the lattice spacing can not matter. However, actually a factor of 3-5 increase was experienced 
"only" , therefore we have to check that the dynamical scaling hypothesis 

C(r, h, t) = X'^iXr, X'^/f'h, X-^H), $(r, h, t) = A^'^l'(Ar, X^^h, X-'H). (13) 

is fulfilled by the solution of Eq.(|l|). The observed dynamical scaling behavior of S,max and h{^max) already suggests 
that our system evolves in the scaling regime. As a direct proof for this, we have rescaled the reduced temperature, 
the magnetic field and the a = dt/dh parameter in such a way that we could expect a factor of 2 increase in the 
correlation length if the first equation of Eq.([T3|) is obeyed. The field $ has been rescaled as dictated by the second 
equation. From the rescaled equations the evolution of the order parameter and of the inverse correlation length has 
been extracted. Before the oscillation sets in the agreement with the expectations based on the scaling hypothesis 
is very convincing. It improves as the system approaches to the critical point. This is a clean argument for the 
independence of our results of the lattice spacing. This conclusion is certainly true until the first maximum of the 
correlation length is reached. For the oscillatory motion a second order dynamics is relevant, therefore the scaling 
behavior based on 2; « 2 should be violated. 

The final question to be discussed in this section is to what extent the observed features of the time dependence 
of ^ depend on the algorithm of its determination. In Fig.^ a typical time evolution is displayed for the three 
characteristic lengths introduced in Section II. Using the rrieff — gi{Ci) relations determined in equilibrium we see 
very good agreement of all three before the system is slowed out of equilibrium. Next all three enter an oscillatory 
regime, with about the same amplitude, but with a certain "phase shift" . 

3.5 



3 
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1.5 
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-0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 0.2 
h 

FIG. 4. Comparison of the time dependence of the correlation lengths determined by three different methods: the "direct" 
correlation length, the average grain size and the inverse mass determined with spectral algorithms. Shown is the evolution 
during the non-thermal crossover passing close to the Ising critical end-point. {L = 32, — 0.15, a = 100, for definitions see 
the text.) 

This observation makes it more difficult to relate the location of the first maxima in h to the freeze-out temperature. 
One has to work with that length which is coupled to the microscopical mechanism producing a certain observable. 
For example, the high frequency part of the cr-field is probably well-represented by a gas of particles with the effective 




7 



frequency determined by the "susceptibility dispersion" method |5[ The decay products of these hard particles will 
reflect this frequency. The softest part of the pion spectra probably comes from the coherent decay of the longest 
wavelength components of the field configuration, therefore a characteristic coherence length closer to the "grain size" 
(method ^) will be followed. 

V. THE EFFECTIVE ORDER PARAMETER DYNAMICS 

In statistical physics noisy first order effective equations are used for the longest wavelength hydrodynamical modes 
to describe relaxation phenomena. These equations have to be written down for the order parameter fields and also 
for composite objects which correspond to densities of conserved quantities From this viewpoint it is not clear 
what could be the foundation for the proposal of BR to write down directly an equation for the relaxation of the 
inverse correlation length. 

We follow the more conventional path and discuss first the effective dynamics of the order parameter. After 
presenting considerable evidence for the validity of our approach to OP, we shall be able to build upon it also a quite 
natural interpretation of the observed behavior of the correlation length. 

It is obvious that one cannot account for the oscillatory behavior of the order parameter experienced after passing 
near the critical end-point just using a first order effective dynamics. This means, that our system actually leaves the 
hydrodynamical regime. Therefore, we propose to complete the effective OP-equation with an "acceleration" term. 
For the effective free energy we use the simplest quadratic form which corresponds to an oscillator potential centered 
at the equilibrium value of the order parameter We emphasize that no noise is introduced. 

The proposed linear equation is of the form 

a-^^{h) + a-^r{h')${h) + m%{h') {^{h) - $eg(/i')) = 0. (14) 

In this fully deterministic equation a "dot" means derivation with respect to h. In view of the relation h = at, the 
derivatives originally refer to the time. The determination of the /i-dependent coefficients nieq = S,'^q , ^eg follows the 
methods described in previous sections. However, their values are taken not at the actual ft., but at a somewhat smaller 
value h' , which corresponds to an earlier equilibrium state. This is the simplest way to incorporate the "slowing out 
of equilibrium" phenomenon into the proposed equation. 

The shift /i — ft' is established by optimizing the agreement with the measured order parameter trajectory. Three 
physical pictures for this shift were tested and benchmarked by the MS deviation of the reconstructed OP trajectory 
from the measured one (<J^). In the first one a global delay parameter is introduced which acts with equal strength 
before and after reaching the neighborhood of the critical point {6^ — 0.0034). In the second version one switches on 
the ft-delay only when ft > (5^ = 0.094). In the third version the delay grows linearly from to a final value which 
was fitted to the data. This latter method produced the smallest MS deviation (5^ = 0.0029), therefore below we shall 
present results obtained with this method. The representative 5"^ values refer to a 32'^ lattice with r, = 0.15, a = 100. 
As expected the fitted slope of the linear shift, ft — ft' = const, x (ft — ftg) clearly increases as r-i gets smaller. The 
value of the constant coefficient changes monotonically from 0.02 to 0.1 while varies from 0.2 to 0.1. 

The reconstruction of the measured OP-trajectory turned out to be rather insensitive to the value of the relaxation 
rate r(ft'). Furthermore, no ft' dependence could be observed. Chosen in the wide range F = 0.01 . . . 0.6, an acceptable 
agreement of the real and the reconstructed OP-trajectory was found. 

A nonperturbative, near equilibrium determination of F was attempted by stopping ft at a certain value hstop and 
fitting the relaxation of $ towards equilibrium to an exponential rate. The estimate for F was found to be in the range 
0.005 . . . 0.025, independently of ri and a. In a model investigated previously we found that during the equilibration 
the relaxation rate approaches its perturbative value strictly from below pl| ] . 

Eq.(^ was solved with the initial values for the OP deduced from the fuU system: l>(ft = -0.2), $(ft = -0.2). OP 
trajectories resulting from different realizations of the initial thermal ensemble give slightly different initial conditions 
for Eq. (p^). This uncertainty sets the error of the reconstructed trajectory. (The average variance of the solution of 
Eq. (|l|) is al^ w 0.0001 for ft < and CTj^ « 0.0006 for ft > 0). 
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FIG. 5. a) Exact evolution of OP and the solution of Eq.(p^) for two values of a. b) Exact evolution of the correlation length 
and its values estimated with help of the effective OP-trajectory (see text) for two values of a. {L = 64, Vi = 0.15) 

In Fig. 1^ we compare the true and the reconstructed OP-trajectories. The quahty of the agreement somewhat 
fluctuates, but its 6^ value is less than 0.02 in the whole L,T,a > 50-region, considered in this paper. In principle 
one could attempt to improve further the analytic interpretation of the OP-dynamics by introducing a memory kernel 
into the equation of $, but we believe that our proposed equation captures the essence of the actual OP-dynamics. 

With this achievement, we can return to the discussion of the non-trivial variation of the correlation length. We 
build our description on our understanding of the order parameter dynamics. 

The evolution of the system is investigated in the close neighborhood of the critical point, in the so called scaling 



13| ) one finds the usual form of the scaling 
ll we have shown the trajectory r = r{h) of 



regime. In equilibrium by choosing the scale factor A = 1/r in Eq. 
behavior of the correlation length and of the order parameter. In Fig. 
the evolution of the system for t — ah. If this trajectory is simple enough, like in the left part of Fig.l, a unique r{£) 
function can be extracted, for instance, from the first relation and after its substitution into the second a functional 
relation i^($) is obtained. One may then deduce the piecewise unique function c?<I>(^)/(i^, which can be used to rewrite 
the effective equation of $ ( p^ ) into an equation for S,{h). 

Motivated by this argument, we have plotted ^(t) against ^(t), both measured during the actual real time evolution 
of the system. It turns out that a unique functional relation ^($) can be recognized in the regime of monotonic 
OP-evolution (cf. Fig. 2a). In the oscillatory OP-rcgime first one has to average the ^ = "^^T// values taken at 
different passes through $. The function ^($) is extracted only after this step. A very good fit valid for the whole 
evolution period of the form 



2.3, 



(15) 



was obtained, with the parameter p slightly depending on the velocity of the /i-variation. The agreement of the 
measured ^(/i) function with the one reconstructed by mapping the computed time evolution of the OP using the 
above relation is quite spectacular (see FigJ|b). 

We conclude this section by discussing a more "theoretical" approach to the determination of F, vrieq and of $eg, by 
observing that these quantities refer to (near) equilibrium situations. In our very simple model their nonperturbative 
values were easily determined from the microscopical data. In case of more realistic models, however, one may attempt 
to use perturbative estimates for the masses, the damping rates and the equilibrium order parameter. We have tested 
the quality of the reconstructed OP-trajectory in the present model also when these coefficients were taken from 
perturbation theory. 

As it is well known even a resummed perturbation theory fails in the vicinity of the critical temperature due to its 
bad behavior in the IR regime. The divergence of the correlation length and as a consequence also of the effective 
expansion parameter AT^ (see for example [Q) excludes its use in the scaling regime. In a finite system, however, 
with an IR cut-off L we could attempt to extract the equilibrium mass, the equilibrium OP value and the damping 
rate of the OP using two-loop lattice perturbation theory. The mass comes from the second derivative of the effective 
potential at the minimum, the OP is the location of this minimum, and for the damping rate we use the formula 
derived in ||l^. The two- loop perturbative effective potential was computed with self-consistent propagators on a 
finite lattice and we extracted from it the position of the minimum and the second derivative in this point along the 
(r, h) route shown in Figj^. The mass used in the propagators was determined from a one-loop gap equation. 

The expression of the two-loop effective potential on lattice reads as 



_ 'J' 



uj{ii, ^) AT2 
f 8L6 



;2(k,^) 
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where ti;^(k,/i) = + /i^. The mass was determined from the one- loop gap equation 

2 2L-^ ^ t^2(k^^) 

In view of Eq. we put A = 6. 

The inaccuracy of the two-loop perturbation theory comes overwhelmingly from the fact that for the lattice spacing 
considered in this paper, the estimate of the critical temperature exceeds its value determined by us numerically at 
mjaa; = 1 by 25%. Repeating both the numerical and the perturbative calculation with Imjaa; = 0.5 the deviation 
diminishes to 5%. Despite all inaccuracies of the perturbation theory the solution of Eq.jT^) based on the perturbative 
potential V2-100P turns out to follow quite closely the real trajectory, though because of the ill-determined OP values 
the deviation has also a systematic error (J^ ss 0.014). 

For an alternative estimate of ^eq{h) we used Widom's scaling form ||l^ like it was done by BR. This yields £,eq{h) 
values close to the measured equilibrium correlation length. The solution of Eq.(p^ fitted to the measured curve with 
this coefficient is only of slightly lower quality than the fully nonperturbatively reconstructed one (5^ « 0.0061). 



VI. CONCLUSIONS 



In this paper we have presented a detailed discussion of the real time non-equilibrium evolution of the classical 
field theory when it passes nearby the critical end point in its (r, h) phase diagram under the influence of a time 
dependent external magnetic field. The numerical investigation was focused on the variation of the correlation length. 
A quick but accurate method for determination of the instant non-equilibrium mass of the <I>-field was employed, 
relying on the parametrical relation of the spatial domain of coherence of field fluctuations to the mass. 

For a wide range of the rate of variation of the magnetic fleld we have experienced a slowing out of the system 
from the hydrodynamical regime. This phenomenon demonstrated itself in the behavior of all quantities used for the 
characterization of the system: the average kinetic energy per site, the order parameter and the correlation length all 
showed oscillations when the value of the external magnetic fleld moved beyond the point of "supercooling" . 

A simple effective equation was shown to describe this dramatic feature of the OP-evolution. The coefficients of 
the second order differential equation written down for the order parameter take their values from equilibrium. They 
can be determined from separate simulations, but also perturbative estimates led to reasonable description of the 
OP-evolution. A more important feature of the effective equation is that the coefficient functions should be computed 
for values of the external magnetic field corresponding to some earlier stage of the evolution. The gradual increase in 
this shift reflects the cumulative effect of the critical slowing down in the proposed effective description. 

The motivation for this investigation came from its possible relevance to the cr-meson dynamics in high energy 
heavy ion collisions ||l|. The weakest point in taking over the above features to QCD is the lack of a quantitatively 
accurate mapping between QCD and the $'*-model in the (T, /x)-plane. Still, our results can be compared in a useful 
way with the scenario put forward in for the evolution near a hypothetical critical QCD end point. 

In the context of heavy ion collisions, it would be interesting to see the effect of a coherently oscillating long 
wavelength cr-background with variable correlation length, on transversal jet quenching. 
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